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Abstract. We calculate the density of states in the half-filled Hubbard model on a Bethe lattice with infinite 
connectivity. Based on our analytical results to second order in t/U, we propose a new 'Fixed- Energy Exact 
Diagonalization' scheme for the numerical study of the Dynamical Mean-Field Theory. Corroborated by 
results from the Random Dispersion Approximation, we find that the gap opens at U c = 4.43 ± 0.05. 
Moreover, the density of states near the gap increases algebraically as a function of frequency with an 
exponent a = 1/2 in the insulating phase. We critically examine other analytical and numerical approaches 
and specify their merits and limitations when applied to the Mott-Hubbard insulator. 

PACS. 71.10Fd Lattice fermion models (Hubbard model, etc.) - 71.27.+a Strongly correlated electron 
systems; heavy fermions - 71.30+h Metal- insulator transitions and other electronic transitions 
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1 Introduction 

The Mott-Hubbard metal-insulator transition is a fasci- 
nating but very difficult problem in condensed matter the- 
ory |1I2I3I4| . One hopes that the basic mechanism can be 
understood with the help of conceptually simple models 
such as the Hubbard Hamiltonian. This model describes 
spin- 1/2 electrons which move on a lattice (band- width 
W = At) and interact locally with strength U. For small 
interactions, U <C W, the model describes a metal with 

\ a finite density of states at the Fermi energy. For large 
interactions, U 3> W, and when there is on average one 

: electron per lattice site (half band-filling), the model de- 
scribes an insulator because it takes an energy of the order 
of A(U ~>W) - 0(W) > to generate a charge ex- 
citation, i.e., the single-particle gap A(U 3> W) is finite, 
irrespective of any symmetry breaking |1I2| . 

These basic insights leave much room for speculations 
on the properties of the transition itself, for example, on 
the precise value of the critical interaction strength U c at 
which the gap opens. Because the Mott-Hubbard transi- 
tion is a zero-temperature quantum phase transition [2], 
it is difficult to access reliably using analytical or nu- 
merical techniques; only exact solutions for all interaction 
strengths can provide conclusive answers. Unfortunately, 
Hubbard-type models can be solved exactly only in one 
dimension |5I6| . and several features in one dimension are 
not generic to three-dimensional systems. 

The phase diagram of a model in three dimensions can 
often be understood from the limit of infinite dimensions. 
For correlated lattice electrons |7I8| . this has led to the 



formulation of effective single-impurity models whose dy- 
namics must be solved self-consistently (Dynamical Mean- 
Field Theory, DMFT) |3I9I1U| . An alternative appr oach, 
the Random Dispersion Approximation (RDA) |2lllj , also 
becomes exact in the limit of infinite dimensions. How- 
ever, an exact solution of the DMFT is not feasible be- 
cause single-impurity models with a Hubbard interaction 
cannot be solved in general. 

Using various kinds of approximations, two different 
scenarios for the Mott-Hubbard transition have been pro- 
posed: 

Discontinuous Transition. 

The gap jumps to a finite value when the density of 
states at the Fermi energy becomes zero at some crit- 
ical interaction strength U c ^'i the gap is preformed 
above U Cl i < U Ct 2, and the co-existing insulating phase 
is higher in energy than the metallic state. 

Continuous Transition. 

The gap opens continuously when the density of states 
at the Fermi energy becomes zero, U c ,i = C4.2 = U c . 

Various approaches to the Mott-Hubbard transition for 
lattices with infinite coordination number yield results in 
favor of one or the other of the conflicting scenarios; for 
a review, see Refs. |2I3| . Refs. m]-^H] give more recent 
treatments. Insightful physical arguments have helped to 
sharpen the physical and mathematical implications of the 
scenario of a preformed gap transition JH1 _ E21 but have 
not been able to resolve the issue. 

All analytical approaches ^7]-[2S] are necessarily ap- 
proximate in nature, and numerical investigations |12|- 
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|28| of the dynamical mean-field equations involve (i) dis- 
cretization, (ii) numerical diagonalization, (iii) interpola- 
tion, (iv) iteration of the self-consistency cycle, and (v) ex- 
trapolation to the thermodynamic limit. The Random Dis- 
persion Approximation requires similar steps in its 
numerical implementation, apart from step (iv). There- 
fore, the region of applicability of the available numerical 
techniques is not a priori clear either. 

In order to address this issue, two of us (E.K. and F.G.) 
have developed a strong-coupling perturbation theory at 
zero temperature |29| which provides a benchmark test for 
the available analytical and numerical techniques for the 
Mott-Hubbard insulator. We assess the quality of the It- 
erated Perturbation Theory (IPT) and the Local Moment 
Approach (LMA) as well as three numerical methods for 
the Hubbard model in infinite dimensions: the Random 
Dispersion Approximation and two Exact Diagonalization 
(ED) schemes for the Dynamical Mean-Field Theory. As 
we shall show, this comparison will provide insight as to 
the limitations of the various methods and will motivate 
an improved ED scheme, the 'Fixed-Energy Exact Diago- 
nalization' (FE-ED). This scheme successfully reproduces 
our findings from the t/U expansion even down to the 
collapse of the Mott-Hubbard insulator. 

In this work, after the introduction of some definitions 
in Sect. |2 we reanalyze our strong-coupling perturbation 
theory to second order in t/U for the Bethe lattice with 
infinite connectivity. In Sect.|3|we estimate that the Mott- 
Hubbard insulator is destroyed at t/ c sc /i = 4.40 ± 0.09. In 
Sect.rjjwe set up the Fixed-Energy Exact Diagonalization 
scheme for the Dynamical Mean-Field Theory. The FE- 
ED very well reproduces our analytical findings, and con- 
firms our conjecture for the critical interaction strength, 
jjFE-ed/£ _ 4.43 ±0.05. In Sect. we favorably compare 
our results with those from the Random Dispersion Ap- 
proximation (RDA) which provides an independent check 
on our results from the strong-coupling expansion and the 
FE-ED. Therefore, we are confident that we accurately 
describe the Mott-Hubbard insulator for all interaction 
strengths. 

In Sect. El we compare our results from Sects. EJ and 0] 
with two analytical approximations for the Mott-Hubbard 
insulator, the Iterated Perturbation Theory (IPT) and the 
Local Moment Approach (LMA). Both of them provide a 
reasonable description down to U — 6t (IPT) and U = 5.5t 
(LMA), respectively, where the LMA is generally superior 
to IPT. However, both approximations underestimate the 
stability of the Mott-Hubbard insulator, t/'^ T = 5.2t and 

In Sect. [7| we critically examine two numerical ap- 
proaches. For finite systems, the Caffarel-Krauth exact- 
diagonalization scheme 14 26 1*3] does not provide unique 
solutions of the self-consistency equations. Therefore, nei- 
ther the size of the gap nor the shape of the density of 
states can be extrapolated reliably. The exact diagonal- 
ization of the effective single-impurity Anderson model 
in 'two-chain geometry' |rSI27l28| fails because the recon- 
struction of the density of states from its moments is 
a numerically too delicate inverse problem. Conclusions, 



Sect. |H1 and appendices on the IPT and LMA algorithms 
in the limit of strong coupling close our work. 



2 Definitions 

In this section, we discuss the basic properties of lattice 
electrons in the limit of infinite dimensions, and we de- 
fine the Hubbard Hamiltonian, the single-particle Green 
function, and the corresponding density of states. 

2.1 Hamilton Operator 

We investigate spin- 1/2 electrons on a lattice whose mo- 
tion is described by the kinetic energy operator 

r = £v£A- (!) 

where c^ a , c\^ a are creation and annihilation operators for 
electrons with spin a =| , J, on site i. The matrix elements 
iy are the electron transfer amplitudes between sites i 
and j, and fy = 0. Since we are interested in the Mott 
insulating phase, we consider exclusively a half-filled band 
where the number of electrons N equals the number of 
lattice sites L. 

For lattices with translational symmetry, iy = t(i — 
j) and the operator for the kinetic energy is diagonal in 
momentum space, 

f = £e(k)c+ CT c k>CT , (2) 

k;<r 

where 

e(k) = i£t(i-j) e - i ( i -j) k . (3) 

The density of states for non-interacting electrons is then 
given by 

p(e) = ^6(e-e(k)), (4) 

k 

The m-th moment of the density of states is defined by 




and e = t(0) = 0. 

In the limit of infinite lattice dimensions and for trans- 
lationally invariant systems without nesting, the Hubbard 
model is characterized by p(e) alone, i.e., higher-order cor- 
relation functions in momentum space factorize |30| . for 
example, 

Pi2(ei, e 2 ) = — £ S(e 2 - e(k + c&))8(ei - e(k + qi)) 

k 

= <*qi,q 2 <5(ei - e 2 Mei) (6) 
+ (1 - Sq u q 2 )p(e 1 )p(e 2 ) ■ 
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This observation is the basis for the Random Dispersion 
Approximation (RDA) which becomes exact in infinite di- 
mensions for paramagnetic systems, i.e., when nesting is 
ignored |2I11| : see Sect.0 

For our explicit calculations we shall later use the semi- 
circular density of states 



^) = Jrf-(wO ■ (M<W, (7) 



nW/2 

1=1 Aujp Q (uj) 

J-W/2 



(8) 



where W — At is the band-width. In the following, we 
take t = 1 as our unit of energy. This density of states 
is realized for non-interacting tight-binding electrons on 
a Bethe lattice of connectivity Z — > oo [31]. Specifically, 
each site is connected to Z neighbors without generat- 
ing closed loops, and the electron transfer is restricted to 
nearest-neighbors: ii.j = —t/yZ when i and j are nearest 
neighbors and zero otherwise. The limit Z — > oo is implic- 
itly understood henceforth. Later on we shall use either 
of the two equivalent viewpoints, RDA or tight-binding 
Bethe lattice, whichever is more convenient for our con- 
siderations. 

The electrons are taken to interact only locally, and 
the Hubbard interaction reads 



1 

TL\ t 

«,T 2 



1 

*U - - 



(9) 



where ni jCr = cf a Ci. a is the local density operator at site i 
for spin a. This leads us to the Hubbard model 32 , 



H = T + UD 



(10) 



The Hamiltonian explicitly exhibits particle-hole symme- 
try, i.e., H is invariant under the particle- hole transforma- 
tion 



where |i| counts the number of nearest-neighbor steps from 
the origin of the Bethe lattice to site i. The chemical po- 
tential [i = then guarantees a half-filled band for any 
temperature 0- 

For later use we also define the operators for the local 
electron density, h- } = n^f + The operators Sf = 

($A+%M/ 2 > s ! = isfA-tf&M*)' and % = 

(nj i — nj ! |)/2 are the three components of the spin- 1/2 
vector operator Sj. The operator for the total spin S = 
Sj commutes with the Hamiltonian. 

2.2 Green Functions 

The time-dependent local single-particle Green function 
at zero temperature is given by 



(12) 



Here T is the time-ordering operator and (. . .) implies the 
average over all ground states with energy Eq, and (taking 
h = 1 henceforth) 



Ci,a(t) = cxp(i_fft)ci !(T exp(— iHt) 



(13) 



is the annihilation operator in the Heisenberg picture. 

In the insulating phase we can readily identify the con- 
tributions from the lower (LHB) and upper (UHB) Hub- 
bard bands to the Fourier transform of the local Green 
function (n = + ), 



G(w) 



dte^G(t) = GlhbM + GuhbM 



GlhbM = J j Y. \ b t + E ) - ir, 

i.er 

Guhb(^) = — Glhb(— u;) . 



(14) 



The last equality follows from the particle-hole symme- 
try, cf. |HJ. Therefore, it is sufficient to evaluate the lo- 
cal Green function for the lower Hubbard band which de- 
scribes the dynamics of a hole inserted into the system. 

The density of states for the lower Hubbard band can 
be obtained from the imaginary part of the Green func- 
tion H14(l for real arguments via |3~3"] 



1 



■DlhbM = -3G L hb(w) , 



lo + H-Eo )c i<a ) , (15) 



with /^Ljjg < to < Mlhb ^ 0. Particle-hole symmetry re- 
sults in £)tjhb( w ) — D L n B (—Lj) so that the single-particle 
gap in the Mott-Hubbard insulator is given by 



A(U) = 2\^+ HB (U)\ >0 



(16) 



We define the (shifted) moments M n (U) of the density 
of states in the lower Hubbard band via 

M n (U)= dw(w + - Dlhb(o;). (17) 

In particular, from (|15|l we find that [33] 

M (U) = 1 , (18) 

(19) 

are two useful sum-rules which we shall employ later. 



3 Strong-coupling expansion 

In this section we first summarize the Kato-Takahashi 
strong-coupling perturbation theory. Next, we prove that 
the entropy density of the Mott-Hubbard insulator is fi- 
nite to all orders in perturbation theory. We then derive 
the density of states to second order in l/U; for further 
details, see Ref. [35]. Finally, we present explicit results 
for physical quantities including the gap and an estimate 
of the critical interaction strength. 
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3.1 Kato-Takahashi Perturbation Theory 

Based on Kato's degenerate perturbation theory |34|. Ta- 
kahashi developed a perturbation expansion in l/U for 
the Hubbard model at zero temperature |35j . For large 
interaction strength U, the set of ground states \ip n ) of H 
in (I10|) can be obtained from states \<fi n ) without double 
occupancy, 



\^n)=P\<Pn) , P \<t>n) = \<Pn) 



(20) 



where Pj projects onto the subspace with j double occu- 
pancies and r is an operator to be determined. 

Eq. (|2(Jfl is readily interpreted. In the large-coupling 
limit, T in (|10(l is considered to be a perturbation to UD. 
Addressing the ground states we therefore start from ei- 
genstates with zero double occupancies \4> n ) into which 
the operator t successively introduces double occupan- 
cies and holes to generate the ground states \\j) n ) of H. 
The operator t reduces all operators to the subspace with 
zero double occupancies. In particular, r + T = Pq so that 
overlap matrix elements obey (VVnlVVi) = (4>m\4 > n}- 

The Schrodinger equation for the ground states 



leads to 



Hf\<l> n ) = E f\cp n ) (21) 



h\4> n ) = E \<}> n ) , h = r+Hr, (22) 



i.e., the eigenvalue does not change under the transforma- 
tion. The creation and annihilation operators are trans- 
formed accordingly, 



cf = r+c+ r 

1,(7 ljCT 1,(T 



h,<j >-> c- ha = P + c- ha r . (23) 



The derivation of the explicit expression for r can be 
found in Refs. |34I35| . where it is shown that 



r = PPn PnPP 



-1/2 



(24) 



Here 



P = P Q - S ri TS r2 ■ ■ ■ TS rm+1 , 



m=l r\-\-...-\-r- m j r i=m 



s° = -P° , 

s r = ( -1 ) r El = izlHsr 

Jjr jr ~ Jjr 



(25) 



The m-th order term in P contains all possible electron 
transfers generated by m applications of the perturba- 
tion T; its contribution is proportional to (1/U) m . The 
square-root factor in (|24|l guarantees the size-consistency 
of the expansion, i.e., it eliminates the 'disconnected' di- 
agrams in a diagrammatic formulation of the theory |35| . 



The square root of an operator is understood in terms of 
its series expansion, i.e., 



(PoPPo) 1/2 =P +Y 



(2m- 1)!! 
' (2m)!! 



Po(P -P)P 



(26) 

The Green function for the lower Hubbard band becomes 



GlhbH = Y ( c t,a [ w + ( h ~ E o) - W 



where (. . .) now implies the average over all states 
with energy Eo. 

A straightforward expansion in l/U, 



(27) 



~a,„ = £ (-u)- m P ^Po 



rn—0 

u 



h=-(L-2N)+J2u- 



m=0 



gives 



h = PoTPq , 

h t = -PoTSfPo , 

k 2 = PoTSfsfPo + hPho 



for the Hamiltonian and 

~<o) _ , 

$1 = \Jt , 



(28) 



(29) 
(30) 
(31) 



(32) 
(33) 



1 



4^ = Ci,aSTST + TSci, a Sf + ici>i + ^1C;, CT (34) 

for the annihilation operators. In the derivation of l|29l) 
(|3~4"Jl we have used the fact that ci, ff in (|2*7|) acts on states 
\4> n ) with no holes and no double occupancies; see also 
Sect. 13.21 We have also used the fact that in the Green 
function for the lower Hubbard band of a half-filled Bethe 
lattice we can express the Hamiltonian h in terms of a 
polynomial in the bare hole- hopping operator ho 



3.2 Ground states 

To leading order in 1/U, the states \<j> n ) are eigenstates 
of the Hamiltonian h\. Because there are no doubly occu- 
pied sites in \<j> n ) and because the system is half filled, all 

sites are singly occupied ('spin-only states'). Therefore, hi 
can be expressed in terms of spin operators. For general 
hopping amplitudes t(i — j) (t(0) = 0) on a regular lattice, 
one finds 



^=£2|i(i-j)| 2 fsi-Sj-J 



ij 



(35) 
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In momentum space this becomes 



hi = Yl J (q) 



Sq ' S_ q — 6 qj 



L 



^ J ( f i) = Z^ e(k)e(k + q) 

k 

/oo poo 
de x ei / de 2 e 2j o q (ei,e 2 ) 
-oo J — OO 



(36) 



Here S q = \/l/LJ2j exp(— iqj)Sj is the Heisenberg spin 
operator in momentum space. 

Since the joint density of states I© factorizes, the Hei- 
senberg model (|3*5f in the absence of any nesting reduces 
to 

^S 2 L 



hx = J(0) 



L 



(37) 



in infinite dimensions. This shows that all global singlets 
are ground states with energy E^p = —J(0)L/(4U) = 
—e 2 L/(2U), see (0- For the Bethe lattice with the semi- 
elliptic density of states (J7J| we have e 2 = 1 so that Eq = 
-L/(2U), 

To first order in 1/U, the degeneracy of the spin-only 
states is not lifted because, in the thermodynamic limit, 
all states with S z — are also global singlets, S = 0. 
More precisely, the ground-state entropy density is given 
by s = ln(2) + C(ln(L)/L) 0. 

This degeneracy is not lifted to all orders in perturba- 
tion theory. We illustrate the argument in the next non- 
vanishing order of the perturbation theory. As was shown 
by Takahashi in the presence of particle-hole symme- 
try, the Hamilton operator for the ground state to second 
order in 1/U vanishes, and the third-order contribution 
can be cast into the form 



h=Y< 2|i(a — b)| 2 |t(b - c)| 2 (S a -S c 



V t(a - b)i(b - c)t{c - d)t(d - a) 

a,b,c,d;b7^d,a7^c 

5 (s b • S c ) (s a ■ S d ) + 5 (s a ■ S b ) (s c • S d ) 

-5 (s a • S c ) (S b ■ S d ) - i (S a ■ S b ) 
-j(Sb-S c ) -i(s c -Sd)- j(s d -S a ) (38) 



■- (s a • s c ) - - ( s b s d 



1 

To 



We have omitted terms where two sites are connected by 
four hopping processes as they give a vanishing contribu- 
tion in the limit Z — > oo. By Fourier transforming, we 
can ensure that the factorization of the correlation func- 
tions as in © hold. Therefore, in infinite dimensions in 
the absence of perfect nesting eq. (f^ reduces to 




The second term describes the (usually ferromagnetic) 
ring exchange whereas the first term favors ground states 
with spin zero. 

Like hi in l|37|l . /13 in (|39|l is solely a function of the 
operator for the total spin. This argument is readily gener- 
alized to higher orders in the perturbation expansion, i.e., 
h is a function of S 2 where the leading term favors total 
spin zero. Therefore, as long as the perturbation series in 
l/U converges, the degeneracy of the ground state is not 
lifted. As a consequence of its huge degeneracy, s = ln(2), 
each lattice site is equally likely to be occupied by an elec- 
tron with spin f or J,, irrespective of the spin on any other 
lattice site. This is a characteristic feature of the Mott- 
Hubbard insulator in the limit of high dimensions. 

On the Bethe lattice there are no closed loops, and the 
ring exchange is absent, e 4 = 2e 2 = 2. Thus, the ground- 
state energy to third order becomes E^ = —L/(2U 3 ). 

3.3 Green function to second order 

As shown in Ref. 29. , to second order in 1/U the Green 
function for the lower Hubbard band can be expressed in 
the form 



1 



G LH b(v) = i^2{ttAho,U) z + g(h ,U) 



(40) 

where z = lo + U — and s(e, U) and g(e, U) are polyno- 
mials in e and 1/U. They are called the shape-correction 
factor and the gap-renormalization factor, respectively. Up 
to second order in 1/U, the general structure of h in |T2pJ) 
fpTTfl and of Ci.cr in JHUhlElJl implies 



s{e,U) = 1 + 



U 



U 2 



, m . gi,2e + gl,0 . ff2,3€ +ff2.ie 
g(e,U) = e+ + (41) 



U 



ir- 



because each electron transfer operator T generates one 
bare hole-hopping operator ho at most. 



3.3.1 Gap-renormalization factor 

We determine the coefficients gij first. To this end, we set 
s(e, U) = 1 in (|40(l and expand this equation in 1/U . We 
equate the result of this expansion with the corresponding 
expression from (|27|l . To first order in 1/U we find 



E 







-2 r 


z - 


-h 





gl,2(M 2 +51,0 



(42) 
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As we need to fix two coefficients only, we further expand 
this expression in 1/z. Keeping the first two non-trivial 
orders we arrive at 



51,2 + .gi,o = j \ d t 



(43) 



and 



25i,2 + 51,0 = T E (K<r 9i,2{hof + 5i,o(>o 



It is not difficult to calculate the expectation values in 
these equations from the definition (|30|) - Recall that (. . .) 
implies the average over all spin-only states. In eq. I|43[) we 
note that an intermediate nearest-neighbor pair of hole 
and double occupancy can be created with probability 
1/2 by hi everywhere but at site i. In the sum over the 
whole lattice, two contributions are therefore missing. Al- 
together, hi generates the contribution — (L — 2)/2 in 
which, together with the first-order contribution L/2 from 
the ground-state energy, results in 31,2 +51,0 = 1- In EH 
the three-site term in hi contributes with probability 1/2. 
Thus, 251,3+51,0 = 1-1/2= 1/2. Therefore, 31,2 = -1/2 
and 31,0 = 3/2, and hf = 31,2^0 + 5i,o =-(&§- 3)/2. 

In second order, we split the contributions arising from 
hi and those from the correction terms. Let C = l/(z+ho) 

and&i = hi-E^-hf, h,f = 32,3^+52,1^0 = hf a +hf b , 
then 



i.tr 



L 



J. J2(£t*Ch 2 Cc ita ) , (45) 

i,tx 

l^^ta^Ch'idc^) . 

i.tr 

(46) 



Expanding in 1/z, the corresponding equations for 
read 



252, 3:a + 52,l:a = ^ E {K,a^' h ^^) > ( 47 ) 
i,(T 

= T E 0%M(M) 3 - 2h Q )h,.) -(48) 



52,3;, 



Note that in l|48[> the operator (ho) 2 — 1 transfers the hole 
at i into its second neighbor shell. After some calculations 
we find 2g 2 , 3 ;a + 32,1,0 = -3/2 and 32,3^ = 3/4 so that 

52,l;o = -3- 

The correction terms were correctly taken into account 
for the Falicov-Kimball model in Ref. but were er- 
roneously omitted for the Hubbard model. The operator 



h'i describes the motion of a hole by two sites whereby a 
spin-flip is induced. In Q46|l the hole hops from site i to 
some site j where the hole induces a spin-flip. This needs 
to be healed after the hole has made an excursion into the 
lattice so that 



\ E (KM^c La ) = -G Q (z)±£ (ctCiKfCh 

i,£T \,(T 

= -lGo(4E(^i, 



4 
3 1 



(49) 



This leads to hf = (3/8)/i , i.e., 52,3:6 = 0, g 2 s-b = 3/8. 
Altogether, 32,3 = 32,3;a + 52,3;b = 3/4 and 32,1 = 

52,i;a +52,i;6 = -21/8, and "hf = 3h [2(h ) 2 - 7]/8. The 
gap-renormalization factor becomes 

, s e 2 -3 3e(2e 2 -7) , 

g(e,tO= e --^+ \ u2 1 (50) 



2U 

through second order in 1/U. 



3.3.2 Shape-correction factor 

In order to determine the coefficients Sjj we first employ 
the sum rules l|18f) and (12) . We use g(e, U) from the pre- 
vious subsection and find through terms of order 1/U 3 , 



(51) 
(52) 



dep (e)s(e,U) , 
^3 =-/ dep (e)s(e,U)g(e,U) , 



because E /L = -C//4-1/(2l7)-1/(2C/ 3 ) + C(C/- 5 ). This 
immediately gives S2,2 = — ^2,0 from (|51|) . From (|52|l we 
obtain su = — 1. 

For the other coefficients we expand (I4U|I and (|27|1 in 
1/U and in 1/z. The shape-correction coefficients to sec- 
ond order follow from the equation 



=4e(<^ (m 2 -i 



S2.2 



~(2) 



~(i) 



(53) 



3 1 

4 + 2 



We hereby correct a minor mistake in Ref. where 
s 2,2 = 5/4 was reported. 

Altogether the shape-correction factor becomes 



/ s e 9(e 2 -l) 

up to and including second order in 1/U. 



(54) 
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3.4 Physical quantities 

3.4.1 Density of states 

The density of states to order (1/U)° is readily obtained. 
It results from the motion of a single hole created at site i 
which moves through the Bethe lattice via ho and returns 
to i. The hole has not disturbed the spin background after 
its return to i, i.e., its motion appears to be free. Therefore, 
ph(w) = Po{u + U/2) 29 37 . Now that we have expressed 
G(u>) in H4Q(I in terms of ho we can easily read off the 
expression for the density of states, 



D 



LHB 



dep (e)s(e, U)S (to + U/2 + g(e, U)) 



(55) 

The second-order expressions for the gap-renormalization 
factor g(e, U) and the shape-correction factor s(e, U) are 
given in JSUJ) and 

To all orders in the 1/U expansion, the density of states 
increases algebraically near the gap, 



D 



UHB 



A(U) 



, (56) 
(57) 



This reflects the fact that po(e) displays a square-root in- 
crease at the band edges. As long as the shape-correction 
factor s(e,U) remains finite for |e| — > 2, the shape of 
the Hubbard bands near the edges remains the same for 
all U > U c . This behavior is qualitatively and quantita- 
tively the same as in the exactly solvable Falicov-Kimball 
model 133 . 



3.4.2 Single-particle gap and width of the Hubbard bands 

As seen from (|55|) the density of states becomes zero at 

U 



ff(T2,£/) = 



(58) 



Using (|16|) we find for the gap 

A(U) = U + 2g(-2, U) = U - 4 - I - + Q (u~ & j 



The Hubbard bands have the width 



(59) 



3 

H^LHB = WW = ptuB - -"LHB = 4 + 7^72 ( 60 ) 



up to second order in 1/U. The 1/U expansion therefore 
provides a very good guess for the support for the density 
of states as a function of frequency in the Mott-Hubbard 
insulator for all U. We shall profit from this result in our 
numerical calculations. 

The gap appears to converge rapidly as a function of 
1/U. Using the expansion parameter x — 2/U, we may 
write lf59"f as 



A{V) = U-l-\ 



x + -x 2 + Xx 3 + 0(x A ) 
4 



(61) 



Preliminary results from our calculations to third order 
indicate that indeed Awl. 

The critical interaction strength is determined from 
A(U) = 0. If we denote by U^ 1 the critical interaction 
strength at which the m-th order expression for A(U) van- 
ishes, we find using A = 1 

C/ c (0) =4, 

C/ c (1) = 4.236 [5.9%] , 
C/ c (2) = 4.313 [1.8%] , 
C/ c (3) w 4.357 [w 1.0%] . (62) 

The numbers in brackets give the percentage change to 
the result from the previous order. Under the assump- 
tions that all higher-order contributions go in the same 
direction, and decay similarly fast, we conjecture that 



W c = 4.40 ± 0.09 



(63) 



which allows for another 3% increase of U c with respect to 
the third-order result. We shall see in Sects. [3] and 0] that 
this estimate is in excellent agreement with our results 
from numerical calculations. 

This analysis may also be carried out for the exactly 
solvable Falicov-Kimball model where 

€k = 2-828 , 



<% = 2.414 [-14.6%], 
f^ K = 2.242 [-7.15%], 



U 



(3) 



c,FK 



Uc 



FK 



2.167 [-3.31%] 
2 [-8.35%] . 



(64) 



It appears to us that for the Falicov-Kimball model the 
convergence of the critical interaction strength is not quite 
as fast as for the Hubbard model. The relative changes 
are about a factor of three smaller for the Hubbard case, 
and correspondingly, we have reason to expect that for 
the Hubbard model the exact result is within 3% of the 
third-order result. 



3.4.3 Self-energy and momentum distribution 

Up to an overall coefficient, the single-particle density of 
states is identical to the imaginary part of the Green func- 
tion l|15(l . The Kramers-Kronig relation provides the real 
part as 



^LHB 



= P / da/D L HB(w') 



1 



1 



(65) 

The local Green function per spin G a (u) — G(cu)/2 and 
the single-particle self-energy S a {uS) = —U a (—ui) are re- 
lated by 



dep (e)- 



1 



w — e — E a {tjj) + it) sgn(w) 
G>-£,H) (66) 
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in infinite dimensions, and 



G°(Y) 



holds for the Bethe lattice |3J so that 
ui - S a (uj) = G a (uj) + 



G a {u) 



(67) 



(68) 



The real and imaginary part of the self-energy are then 
readily obtained as 



1 



{^G a {oj)f + (3G CT ( W )) ! 



1 



(3£G ff H)2 + (3G CT M) 2 



(69) 



(70) 



Note that this self-energy is not perturbatively linked to 
the weak-coupling limit. For example, its imaginary part 
displays a peak 5(u>) and its real part diverges proportional 
to (— l/o>) as to — > 0. 

The single-particle spectral function A a {t\ lo) is defined 

by 

Aa(e;uj) = -i B gn(w)9f ( ^jrj^) (71) 

_ _i . v gjWM 

^sgn lW j ^ _ £ _ ^^^2 + (CJ^( W ))2 ' 

It does not display a quasi-particle contribution. The mo- 
mentum distribution 



n CT (e) 



^LHB 



(72) 



depends on momentum implicitly via e = e(k). In the in- 
sulating phase, the momentum distribution is an analytic 
function of e. Since n„(e) = 1 — n a {— e) due to particle- 
hole symmetry, the Taylor expansion around e = has 
the form 

n(e) = i - c l£ - c 3 e 3 .... (73) 

For the Bethe lattice we obtain ci = l/(2[/)+e>(E/~ 3 ) 
and, in general, C2„_i = C([/ 1_2 ™). 



4 Exact diagonalization with fixed energy 
interval (FE ED) 

In this section, we first discuss the single-impurity model 
onto which the Hubbard model can be mapped in the limit 
of infinite dimensions. Then we propose the Fixed-Energy 
Exact Diagonalization (FE-ED) as a new scheme for the 
numerical solution of the DMFT equations. The results of 
this approach corroborate our findings from Sect. |3| 



4.1 Dynamical Mean-Field Theory (DMFT) 

In the limit of infinite dimensions [7] and under the con- 
ditions of translational invariance and convergence of per- 
turbation theory in strong and weak coupling, lattice mod- 
els for correlated electrons can be mapped onto single- 
impurity models which need to be solved self-consistently 
|9I10I3| . Unfortunately, these impurity models cannot be 
solved analytically. 

For an approximate numerical treatment, various dif- 
ferent implementations are conceivable; see, e.g., Ref. |31?] 
for a recent implementation. One realization is the single- 
impurity Anderson model in 'star geometry', 



n.-l 



-ffsiAM = 51 e$t;t$°;£ + U [ d i d T ~ 2 ) \ a l a i ■> 



e=l:a 



+ E E ViU+A+d+^cr; 



(74) 



where Ve are real, positive hybridization matrix elements. 
The model describes the hybridization of an impurity site 
with Hubbard interaction to n s — 1 bath sites without in- 
teraction at energies tt with < ei < 62 < . . . < £(n s -i)/2- 
In order to ensure particle-hole symmetry, we have to set 
U = —e n ,-l and V e = V na -l for 1= (n s + 1)/2, . . . , n s - 1. 
Moreover, since we are interested in the Mott-Hubbard 
insulator, we only use odd n s so that there is no bath 
state at e = 0. 

For a given set of (n s — 1) parameters (ee,Ve) the 
model (I74|l defines a many-body problem for which the 
single-particle Green function 



G<?-\t) = -i(f \i(t)dt 



SIAM 



(75) 



can be calculated numerically for n s < 15 with the (dy- 
namical) Lanczos technique. In l|75(l (. . .)siam implies the 
ground-state expectation value within the single-impurity 
model. Typically, the imaginary part of the Green function 
displays n s peaks with large weight and many other small 
peaks whose number depends on the number of states til 
kept in the Lanczos diagonalization. 

Ultimately, we are interested in the limit n s — > 00 
where the hybridization function 



n.-l 

E 

1=1 



v? 



lu — ee + iTysgn(w) 



(76) 



should smoothly approach the hybridization function of 
the continuous problem, 



H(oj) = lim H (ns) {uj) 



(77) 



Correspondingly, the Green function should fulfill 

G<j(lu) = lim G£-\u). (78) 
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At self-consistency, the Green function of the impurity 
problem describes the Hubbard model in infinite dimen- 
sions. As shown in on the Bethe lattice the hybridiza- 
tion function must obey the simple relation 



1 < I < (n s — l)/2, we can be sure that our resolution 
of the Hubbard bands becomes increasingly better as n s 
increases. 

We address the second problem in the following way. 
We apply a constant broadening of width SW to the indi- 

(79) 

vidual peaks in QGa 1 ^ (uj r ), r — 1, . . . , n^. Then, we col- 



This equation closes the self-consistency cycle for the con- 
tinuous problem: we have to choose bath energies and hy- 
bridizations in such a way that the single-particle Green 
function and the hybridization function fulfill (|79|) . 



4.2 Implementation of the FE-ED 

Any numerical scheme for the DMFT faces two problems. 
First, and foremost, the DMFT requires the Green func- 
tion G g (uj) on the whole, continuous frequency axis but 
only the information for (n s — 1) energies eg is provided. 
This is in contrast to standard numerical problems in 
many-body theory where a given energy interval, typically 
of the order of the band-width W, is resolved to accuracy 
W/n s . Concomitantly, in numerical DMFT calculations it 
is not a priori clear how observables scale as a function of 
l/n s . Moreover, there can be more than one self-consistent 
set of parameters {eg, Vg) for fixed n s . 

Second, the self-consistency condition l|79l) holds for 
the continuous H(oj) and G a (uj), not for their discretized 
counterparts. Therefore, the various ED methods will dif- 
fer in the way the information is extracted from G^^iuj) 
in one iteration in order to specify the input parameters 
(eg,Vg) for the next iteration. This is another source of 
ambiguity. Since the scaling as a function of l/n s is not 
clear, it cannot be guaranteed that different schemes will 
ultimately coincide as n s — * oo. 

In this work, we use the results of the 1 /U expansion to 
circumvent the first problem. We have seen in Sect. El that 
the Hubbard bands are distributed symmetrically around 
uj = with width W LH b = TU UH b = W* « W. Thus, 
we actually need to resolve a finite frequency interval. In 
practice, we use W* — 4.5 as our maximum band-width. 
Then, we determine the onset of the upper Hubbard band 
self-consistently, A(U)/2; see below. To this end, we start 
with some input guess A- lri {U). In practice, we use the 
second-order estimate of the gap from l(5§)l to speed up 
the convergence of this outer self-consistency cycle. 

We choose to discretize the Hubbard bands equidis- 
tantly, i.e., we fix the energies eg for 1 < I < (n s — l)/2 



by 



A in (U) 



SW 



SW 



2W* 



1 



(80) 



For n s < 15 and U > W none of the eg can move outside 
the Hubbard bands. By fixing the energies at the centers 
of the equidistant intervals 



(£- 1)SW, 



lect the weight into the intervals Ig and assign this weight 
to wg = Vg . In the Caffarel-Krauth scheme |2fill8| a ficti- 
tious temperature and a x 2 -fitting procedure is used. We 
have found that our simple approach is equally suitable. 

Typically, the weight of the peaks at energies outside 
the Hubbard bands is very small. Thus we set 

dLuJ2^G^(LU r ) 

-I r=l 

0{uj - uj t + SW/2) - 0{u - oj r - SW/2) 



SW 



.(82) 



esw] 



(81) 



In order to generate an educated input guess for Vg — ^fwg 
we use the results from second-order perturbation theory 
in PI). 

Now that we have determined the input energies eg and 
hybridizations Vg, we can start the inner self-consistency 
cycle. We perform a dynamical Lanczos diagonalization 
which gives QG^^^lo,,) for Lanczos frequencies. Most 
of the peaks carry little weight, only n s of them are sig- 
nificant. Using (|82|l . the Green function provides the hy- 
bridization elements for the next iteration of the inner 
self-consistency cycle. We have checked that, for fixed n s , a 
unique solution for wi is found for various starting choices. 

After we have determined the parameters Vg self-con- 
sistently for fixed n s and a given U , we may calculate the 
single-particle gap from the difference in the ground-state 
energy at half band-filling, Eq 1b \N; U), and the ground- 
state energy for one particle more than half band-filling, 
E ( ns) (N+ l;U), as 

A E > (n °\U) =2(E( ns) (N + l;U)~E { Q ns) (N;U)) . (83) 

After the extrapolation to the thermodynamic limit, 

A E {U) = lim A E ^(U) , (84) 

n 3 — >oo 

we obtain a new estimate for the onset of the Hubbard 
bands, and we could restart the outer self-consistency cy- 
cle. Unfortunately, the extrapolation n s — > oo has limited 
accuracy, and it is difficult to obtain a truly converged so- 
lution which fulfills A E {U) — Ai n (U) ('energy criterion'), 
even at very large values, e.g., U = 14. An example for 
U = 6 is shown in Fig.^ It is seen that the results around 
A[ n — 0.88 are almost but not quite converged. 

In order to stabilize the outer self-consistency cycle, 
we need an additional criterion. To this end, we confirm 
numerically that the density of states near the gap in- 
creases algebraically; see l|5l)| . The weight w\ = V? of the 
peak in Dtjhb^) at e\ represents the integral of the den- 
sity of states from A(U)/2 up to e x = A(U)/2 + 0(l/n s ). 
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0.05 0.1 

II n s 

Fig. 1. Gap A E - (ns) (U) from for U = 6 as a function of 
1/ri.jj for various A- m . 



8 1 1.2 1.4 1.6 1.8 

Fig. 2. Weight of the first peak in the upper Hubbard band 
for n s = 7, 9, 11, 13, U = 6, and various Z\m. For a square-root 



increase of the density of state, 



2/3 



should scale linearly with 



l/n s , see 1851 . The intersection of the extrapolated curves give 
A w ( 6) . Inset : expansion of the region around the zero intercept . 



Using (|5oT) we should find 



A W (U) 



(85) 



The linear behavior with cvfe-ed = 1/2 is confirmed in 
Fig. [3 for [7 = 6. The plot also provides A W (U) for a given 
A- ln {U) as the intersection of the extrapolated curves with 
the frequency axis ('weight criterion'). 

The two values for the gap from the 'energy crite- 
rion' l|84(l and from the 'weight criterion' (|85[) do not agree 
perfectly. Therefore, we use the self-consistency condition 

^i^. ^B , (86, 

which determines A FE ~ EI) (U) after convergence of the 
outer self-consistency cycle. The difference, 



s 



1- 



o 



T 



FE-ED 
Second order 



6A 



FE-ED 



(U) 



\A E {U)- A W (U)\ 




(87) 



is our estimate for the accuracy of Z\ FE_ED ({7). We call 
this implementation of the DMFT the 'Fixed-Energy Ex- 
act Diagonalization (FE-ED)' scheme. 

The converged FE-ED gap as a function of the inter- 
action strength is shown in Fig. |3J The agreement with 
the 1/U expansion is excellent for U > 5 because the 
two curves agree within the error bars. Unfortunately, the 
uncertainty in SA FE ^ ED (U) increases towards the tran- 
sition. Using a cubic polynomial extrapolation for data 
points U > 4.8, as shown in the figure, leads to U F E ~ ED = 
4.44. As seen from the figure, the spline interpolation also 
touches the data point for U — 4.6, to within its error 



Z. 



Fig. 3. Mott-Hubbard gap as a function of U in FE-ED in 
comparison with the 1/U expansion to second order 15911 . The 
error bars give the two values which enter 1861 as obtained 
from the energy and weight criterion. The dashed line is a 
cubic polynomial interpolation based on the data points for 
U > 4.8. 



bars. If this data point were included, the critical value 
would be U FE ~ Erj = 4.42. Thus, we give 



jtFE-ed = AM ± aQ5 



as our best estimate for the location of the closing of the 
gap. 
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CO 

Fig. 4. Density of states, D a (uj) = D(u>)/2, as obtained from 
the 1/(7 expansion, eq. I|55|l . and the FE-ED for n 3 = 15; (a) 
U = 6, (b) U = 5. 



In Fig. 01 we compare the density of states of our FE- 
ED for n s = 15 with those of the 1/U expansion to second 
order. Both for [7 = 5 and U = 6 the agreement is very 
good. Our FE-ED thus allows us to determine the density 
of states quantitatively with increasing resolution as n s 
increases. 

Aside from the fact that the 1/U expansion provides 
an educated guess for the position of the Hubbard bands, 
the FE-ED is an independent numerical method based on 
the DMFT. The good agreement of the results mutually 
supports the applicability of either method. Therefore, we 
are confident that neither of the limitations poses a se- 
rious problem for our investigation of the Mott-Hubbard 
insulator, i.e., the problem of convergence of the series 
expansions in 1/U and in l/n s , respectively, are not too 
serious in our approaches. 



5 Random Dispersion Approximation 

In this section, we present numerical results from the Ran- 
dom Dispersion Approximation (RDA). This method is 
not based on the DMFT and thereby provides yet an- 
other independent check for the validity of our results as 
obtained in Sects. 01 and 0] 

In the Random Dispersion Approximation, the disper- 
sion relation e(k) in the kinetic energy is replaced by a 
random quantity e RDA (k) where the bare density of states 
acts as probability distribution, 



^) = |E^- eRDA ( k ))- 



(89) 



This is the characteristic property of the dispersion re- 
lation in infinite dimensions [2] so that the RDA with 
the semi-elliptic density of states (JTJ becomes exact for 
the Bcthe lattice with infinite coordination number in the 
thermodynamic limit L — > oo. 

All correlation functions factorize, analogously to 
For the mth-order correlation function, 



Pi...m(ei 



1 m 

lEIP ^-e(k + qO) , (90) 



k l=l 



we give a recursive algorithm. Let pi(ei) = p{^i) be the 
bare density of states, for m = 2, see ©, and for m > 3 
we define 

m 

Pl...(n-l);n,...,m( e l' • • • J £ ™) = Y\. ^ _ ^liilj) 

i,j=n 

m n— 1 

nn^-w ( 9i ) 

i—n j — 1 

Pl...m( e l) ■ ■ ■ i e m) ■ 

In this quantity, all q; for I > n are pairwise different, and 
they are also different from those with i < n — 1. The 
recursion formula 

Pl...n;(n+l),...,rn( e lj ■ ■ ■ J e m) = 

(ei, . . . , e m ) 

n-1 

1 - II ( X - 5 ^{e n - e,)) (92) 



i=i 



<Pl...(n-l);(n+l),...,m( e l 



i e m) 



starts with /fi... m (cii • ■ ■ Em) = Pi...m-. (ei, . . . e m ) and ter- 
minates at 



Pl;2,.,m(ei, •••em) = n^ ej ) II ( X 6 1i,1j] 

1=1 i,j=l 

For example, 



(93) 



Pi23(ei,£2,£3) = 1 J]^ 1 ~~ S i3,n S ( e 3 - £()) Pi2(ei,£2) 



1=1 



+Pi2;3(ei,e2,£3) 



(94) 



with 



Pi2;3(ei, £2, £3) = Sq u q 2 S{ei - e 2 )p(ei, e 3 ) 

+P(ei)p(£2)p(e 3 ) (95) 

X (! - 5 qi,q 2 )( 1 - ^U.qaX 1 ~ <*q2,qs) ■ 

In order to put the RDA into practice, we choose a 
one-dimensional lattice of L sites in momentum space 
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and determine the dispersion relation e(fc) as the solution 
of the implicit equation 

k/2 = (2e{k)/W) [1- {2e(k)/W) 2 } 1/2 + &rcsm{2e(k)/W) . 

(97) 

This choice guarantees p(e) = po(e) in the thermodynamic 
limit. 

Next, we choose a permutation Q a for each spin di- 
rection er which permutes the sequence {1,...,L} into 
{Q<t[1], ■ ■ ■ , Q<t[£]}; in this way there are (L\) 2 indepen- 
dent realizations for a given finite L. A specific realization 
of the RDA dispersion is then defined by Q = [Qj, QJ. 
Then, the numerical task is the Lanczos diagonalization 
of the Hamiltonian 



(98) 



6=1 



In this way we obtain the ground-state energy E^(N; U), 
the single-particle gap 

A Q (U) = 2[E^(N = L + 1; U) - E®(N = L; U)] , (99) 

and the momentum distribution 

> s (e;C/) = 5X;<*W> fuy , (ioo) 



n - 



where (. . .) denotes the ground-state expectation value for 
the realization Q. 



5 
4 
3 
2 
1 




g 
< 



1 1 1 1 II 1 1 1 II 1 1 1 II 1 

- o t/=3.82 

- ■ £7=5.09 

A U=631 ^ 


1 1 1 1 1 1 1 1 1 










i 1 i i i i 1 i i i i 1 i i 


■ i 1 i i i i 1 i 


) 0.05 0.1 0.15 


0.2 0.25 



1/L 



Fig. 5. Gap as a function of system size (L < 15) in the 
Random Dispersion Approximation for various values of the 
interaction strength. The extrapolations are done for even and 
odd system sizes separately. 



As a next step, we obtain all physical quantities for 
fixed system size L by averaging over the number Nq of 
realizations Q. Typically, we choose at least Nq = 100 
for 6 < L < 16. For the physical quantities we obtain 
gaussian-shaped distributions for which we can determine 
the average values, e.g., 



Nq q 
Nq Q 

with accuracy 0(1/Nq). 

3l ' ' ' ' 



(101) 
(102) 



2- 



g 



1- 



0, 



— •— RDA, L even 
RDA, L odd 
Second order 




J I I L 



3.5 4 4.5 5 5.5 6 6.5 7 
U 

Fig. 6. Gap in RDA as a function of the interaction strength. 
Filled symbols: extrapolation from even L, open symbols: ex- 
trapolation from odd L. Also shown is the result from second- 
order perturbation theory in 1/(7 I59H . 



In order to improve slightly the quality of our distribu- 
tions, we impose a filter on our randomly chosen permu- 
tations. For a truly random dispersion, L|£ RDA (i?)| 2 = e 2 
is independent of £; see IpIrJl) . Therefore, we discard those 
realizations for which 



L-l 



> d L 



(103) 



with, d,L ~ 0.2 for L < 16. In this way we admit about 
every second of the randomly chosen configurations. Note 
that there are of the order of (LI) 2 different realizations 
so that our filter does not introduce any unwanted bias. 

In Fig. we show the gap as calculated from l|9"9"|) . 
Typically, the statistical errors are of the symbol size or 
smaller. Due to the sizable odd-even effect, we extrapolate 
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the data for odd and even L separately. This gives another 
check on the influence of finite-size effects. It is seen that 
the behavior as a function of 1/ L is quite regular, i.e., the 
finite-size gap behaves as in generic many-body problems. 

In Fig. we show the result of an extrapolation of the 
finite-size data for even and odd system sizes. We consider 
the agreement with the l/U expansion as quite acceptable 
for U > U c which we can estimate as 



U, 



RDA 



4.0 ±0.4 . 



Finite-size effects result in a smearing of the gap at the 
transition, and a gap exponent cannot be determined. 
However, the fact that the extrapolated gap turns slightly 
negative can be taken as an indication that the gap opens 
at U c with an exponent not too far from unity. 



1 
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Fig. 7. Momentum distribution in RDA for U = 16/ir « 5.1 
for the two largest system sizes. The full line is the result from 
first-order perturbation theory in l/U 1731 . 



Lastly, in Fig. we show the momentum dispersion 
within the RDA for U — 16/ir as 5.1, together with the 
result of the l/U expansion. Obviously, there are large 
finite-size effects which make a finite-size extrapolation 
difficult. The reason for this behavior is related to the fact 
that all finite-size systems appear to have a finite disconti- 
nuity at the Fermi energy as well as a finite single-particle 
gap. For insulators, the apparent jump in n(e) only van- 
ishes in the thermodynamic limit . The agreement be- 
tween the RDA data and the momentum distribution as 
obtained from the l/U expansion is good enough to argue 
that the RDA will eventually converge to the l/U result 
for L — > oo. 

The RDA becomes exact for lattice electrons in high 
dimensions. In contrast to the FE-ED it is not based on 



the DMFT self-consistency equations. Moreover, it does 
not require the convergence of the l/U expansion. Despite 
the observed limitations in accuracy, the results from RDA 
confirm of our findings in Sects. |21 and 0] the l/U expan- 
sion and the FE-ED provide an accurate description of the 
Mott-Hubbard insulator whose gap opens at the critical 
interaction strength U c = 4.43 ± 0.05. 



(104) 5 Comparison with analytic approximations 



In this section we compare our results with two analytic 
approximations to the DMFT which become exact for the 
Bethe lattice in the strong-coupling limit, i.e., they re- 
produce -Dlhb(w) to order (1/U)°. The Local Moment 
Approach (LMA) outlined in appendix IB1 fulfills this 
condition by construction, whereas this is a non-trivial re- 
sult within the Iterated Perturbation Theory (IPT) |4U|: 
a proof |41j is given in appendix lAl 

The Hubbard-III approximation 02 a ^ so fulfills this 
criterion but it fails to reproduce qualitatively the density 
of states for smaller interaction strengths [2H|. Therefore, 
we do not further discuss the Hubbard-III approximation. 



— 6.1 Iterated Perturbation Theory 



For the Z — > oo Bethe lattice, the self-energy and the local 
single-particle Green function are related by 



G a {u) = 



1 



uj - G<j{uj) - ' 
In IPT, the self-energy is approximated by 

r°° d<? 

s (7 (oj) = u 2 ^±n{Q)g a (oj-Q) 

J-oo 27r l 

Here Q a {ui) is the host Green function, 

1 



UJ - G a {uj) 1 

and we have defined the polarization bubble 

n(uj) = -( ^rG<r{U\)G<T{U\ - w) 
J-oo 2tti 



(105) 



(106) 



(107) 



(108) 



in terms of the host Green function. 

It is helpful to consider the low-frequency behavior of 
the IPT equations. First, consider the host Green func- 
tion (|107fl . whose low- frequency behavior in the insulator 
is given generally by 



Go (w) 



too 



uj -> 



uj + ir] sgn(w) 
where the pole weight |mo| is given in terms of 

dG a (uj) 



I7i I 



doj 



(109) 



(110) 
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as 



KHa + Nr 1 . (in) 

Q a (oj) also contains a Hubbard-band contribution, which 
is strictly separated from the pole at oj — throughout 
the insulating phase, whence we may write 



\m \ 



(112) 



oj + ii]sgn(oj) 

which actually defines Q^ and (oj). In IPT the low- frequency 
behavior of Q a (oj) directly determines the low-frequency 
behavior of S„{u); inserting 11121) into (|106[> and (|108|l 
yields 

S a (oj) = - oj^O, (113) 



where 



TT2 

A = ^|m | 3 



To satisfy self-consistency (|1L)5[) requires 

1 



.4 



171 1 



(114) 



(H5) 



so that (|lllfl , l|114|) and i|115fl determine the low-frequency 
behavior of the problem. In particular, they yield 

(116) 



M a/ 1 l m d 



= u . 



The denominator 
value when 



m 



mo 



reaches its maximum 



\m \ 



(117) 



thus locating the destruction of the Mott-Hubbard insu- 
lator at 



U, 



IPT 
c.l 



= 3^3 re 5.2 , 



(118) 

as earlier reported by Rozenberg et al. [21! • The IPT pre- 
diction lies considerably above our best estimate U c = 
4.43 ±0.05. 

This discrepancy can already be seen for larger inter- 
action strength where 



A^(U)=U-A-^- 



(119) 



The first-order correction can be obtained analytically, as 
shown in appendix^ It is more than twice as large as our 
exact coefficient. Therefore, the agreement is good only for 
U > 6. The gap in IPT is compared to our FE-ED and 
second-order results in Fig. EI As seen from the figure, and 
shown analytically in |41j , the gap closes continuously but 
with an exponent of 7ipt = 1/2. 

Despite the deviations in the gap values, the IPT cor- 
rectly reproduces the overall shape of the density of states 
for U > 6, apart from the artificial high-energy tails of the 
IPT Hubbard bands. The threshold exponent for the den- 
sity of states is correct, aiPT = 1/2 in the Mott-Hubbard 
insulator. Nevertheless, the deviations from the second- 
order result are quite noticeable at U — 6 which is 
shown in Fig. |5| We conclude that the IPT provides a 
quantitatively correct description of the Mott-Hubbard 
insulator down to U ps 6. However, IPT seriously under- 
estimates the stability of the Mott-Hubbard insulator. 



s 

< 



1- 



Second order 
FE-ED 
LMA 
IPT 








Fig. 8. Gap as a function of interaction strength from second- 
order perturbation theory in 1/(7 15911 . the Fixed- Energy Ex- 
act Diagonalization (FE-ED), the Local Moment Approach 
(LMA), and Iterated Perturbation Theory (IPT). 
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Fig. 9. Density of states for U — 6 in Iterated Perturbation 
Theory (solid line) and in second-order perturbation theory 
in 1/U l!55t (dashed line). 



6.2 Local Moment Approach 

Details of the LMA are given in Ref. As shown in 
appendix iBl the LMA gap to first order in 1/U is given 

by 

3 



A 



LMA 



(120) 



The first-order correction is close to the exact one, and, 
correspondingly, the agreement with our strong-coupling 
result (|59|l is very good down to U ps 5.5. For smaller in- 
teraction strengths, the deviations become noticeable and 
the LMA predicts the collapse of the insulator to occur at 
jjLma _ 4_g2 [53] which is larger than our best estimate 
U c = 4.43 ± 0.05. The LMA yields a critical gap exponent 
7lma = 1 |41j . in contrast to IPT where 7ipt = 1/2. The 
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gap in LMA and our second-order and FE-ED results are 
shown in Fig. |H1 

It is seen that the LMA gap is more accurate than the 
gap from IPT. In particular, this holds true for the density 
of states as the LMA very well reproduces our second- 
order result (|55(l . The comparison for U = 4V2 w 5.66 
is shown in Fig. ^] The deviations are largest around 
the gap. The overall agreement, however, is excellent. In 
particular, the threshold exponent is qlma = 1/2. 



t 1 1 1 1 1 1 1 1 1 1 1 1 1 r 




CO 

Fig. 10. Density of states for U = 4\/2 w 5.66 from the Local 
Moment Approach (full line) and second-order perturbation 
theory in 1/U l|55jl (dashed line). On the scale of the figure, the 
gap from FE-ED and from second-order perturbation theory 
are almost the same. 

The LMA thus provides a qualitatively correct descrip- 
tion of the Mott-Hubbard insulator. It is a quantitatively 
reliable approximation down to U « 5.5. 



7 Other numerical approximations to the 
DMFT 

In this section we discuss two other exact-diagonalization 
approaches. In both of them a discretized version of the 
single- impurity Anderson model (SIAM), Sect. 14. ll is in- 
vestigated in which the site energies ei are determined self- 
consistently. This is the decisive difference to our 'fixed- 
energy' algorithm as presented in Sect.^ 

In this work, we do not discuss the Numerical Renor- 
malization Group (NRG) approach to the Mott-Hubbard 
insulator. The NRG is custom-tailored for the description 
of the metallic phase in the sense that it tries to resolve 
structures around uj = whereas features at high ener- 
gies are broadened on a logarithmic scale ^5]. Therefore, 
the Mott-Hubbard insulator does not display a clear gap. 
Consequently, as can be seen from Fig. 2 in Ref. the 
density of states is finite in the gap, and the NRG Hub- 
bard bands display sizable high-energy tails. Therefore, 



the gap and the width of the Hubbard bands as a func- 
tion of the interaction strength cannot be easily deduced 
from this approach. 



7.1 Exact diagonalization in 'two-chain geometry' 

In contrast to Sect. 14. ll an alternative formulation of the 
Hamilton operator for the SIAM describes two chains of 
length (n a — l)/2 which represent the upper and lower 
Hubbard bands. They are coupled to the impurity site at 
the respective chain origins |3I27| 



#SIAM-tc — 




(121) 



(n 3 -3)/2 

Here particle-hole symmetry results in q > for £ = 
l,...,(n s - l)/2, and V t > for £ = 1, ... , {n s - 3)/2. 

The main advantage of the two-chain geometry lies in 
the fact that on the Bethe lattice the self-consistency con- 
dition {25 

can be implemented exactly even for finite n s . 
In the continued-fraction expansion of the Green function 
for the Hubbard bands l|14fl . the local energies q and the 
electron transfer amplitudes Vg appear as diagonal and 
off-diagonal coefficients, 

1/2 

Glhb^M = — „,„ , T . , 2 /r t ■ (122) 

uj + (7/2 - {Viy/[uj - ei - . . .\ 

The (dynamical) Lanczos procedure directly provides the 
continued-fraction coefficients so that the output of the 
exact diagonalization procedure yields the input for the 
next iteration of the self-consistency cycle without a col- 
lection of weights as, e.g., in Q82J1. 

For this reason, the two-chain ED is rather appealing. 
Moreover, using the results for n s = 3,5,7 and some un- 
specified extrapolation in l/n s , the gap in was found to 
be almost linear as a function of the interaction strength 
above U c = 4.30, in surprisingly good agreement with our 
best estimate (JHHJ - The corresponding results for n s = 
3, 5, 7 are shown in Fig. ^2 If we take only those three data 
points into account and assume a linear scaling in 1 /n s , we 
obtain the gap as shown in Fig. ll2l which is very similar to 
the result reported in [3], and is in reasonable agreement 
with the results from the l/U expansion. Thus, one would 
be tempted to conclude that the method is quite reliable. 
Unfortunately, such a conclusion is not warranted. 

First, bigger system sizes need to be analyzed because 
n s — 3,5,7 correspond to chain lengths of L c = 1,2,3. 
When we follow the original algorithm as described in |2Zj , 
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Fig. 11. Gap in exact diagonalization of the single- impurity 
model in two-chain geometry for n s —3,5,7 for various inter- 
actions U. The lines are a linear extrapolation to l/n s — > 0. 




u 

Fig. 12. Extrapolated gap (see Fig. llH in the exact diagonal- 
ization of the single-impurity model in two-chain geometry as 
a function of the interaction strength, as compared to second- 
order perturbation theory in 1/U 1591 . 



Fig. 13. Gap in exact diagonalization of the single-impurity 
model in two-chain geometry for n s < 13 for various interac- 
tions U. The lines are a cubic spline extrapolation to l/n s — ► 0. 
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Fig. 14. Extrapolated gap (see Fig. I13B in the exact diagonal- 
ization of the single-impurity model in two-chain geometry as 
a function of the interaction strength, as compared to second- 
order perturbation theory in 1/U 1591 . 



we find that the insulating phase becomes unstable at 
U = 8 for n s = 13 (U = 6.4,5.0,3.5 for n s = 11,9,7). 
Obviously, this is a numerical artifact which is not related 
to the metal- insulator transition. We have found that this 
problem can be cured by incorporating the particle-hole 
symmetry into £fsiAM-tc <|121|1 right from the start; this 



was done only 'on average' in Ref. [27|. In the following, we 
show results obtained using this particle-hole symmetry. 

When we include bigger system sizes, we obtain a be- 
havior quite different from the one seen in Fig. 1121 The 
result for the gap as a function of l/n s for n < 13 is 
shown in Fig. El The behavior of the gap as a function of 
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l/n s can no longer be approximated by a linear fit. The 
gap is not regular as a function of 1 /n s but seems to satu- 
rate. When we interpolate using a cubic spline, we obtain 
the gap as a function of the interaction strength as shown 
in Fig. ^] Obviously, the results are rather poor as the 
large-f/ limit cannot be reproduced properly and there is 
nowhere agreement with the results from the l/U expan- 
sion. It is thus seen that the agreement between the data 
reported in |3I28| and ours is fortuitous. 




CO 

Fig. 15. Density of states as obtained from the two-chain exact 
diagonalization for (7 = 6 and n s — 13 (full line), as compared 
to second-order perturbation theory in l/U 15511 (dashed line). 



The poor quality of the two-chain ED is also seen in 
the density of states, as shown in Fig.EDf° r n s = 13 and 
{7 = 6. For all n s > 7 we obtain a four-peak structure 
which does not allow the recovery of the true density of 
states. 

The reason for the failure of the method lies in the 
setup of the self-consistency scheme. In the continued- 
fraction expansion one optimizes the moments of the den- 
sity of states. However, the reconstruction of the density of 
states from its moments is a numerically very delicate in- 
verse problem. In our case, the moments of the upper Hub- 
bard band can be approximated very well by three main 
peaks and many small peaks at (much) higher energy. For 
this reason, the structure of the density of states is essen- 
tially unchanged as a function of n s . The convergence as 
a function of system size saturates, and a reliable extrap- 
olation scheme is not evident. We thus conclude that the 
two-chain ED is not a suitable tool for the investigation 
of the Mott-Hubbard insulator. 



7.2 Cafrarel-Krauth exact diagonalization in 'star 
geometry' 

As an alternative to the two-chain discretization of the 
previous section, Caffarel and Krauth used the discretiza- 
tion of the single-impurity Anderson model in 'star geom- 



etry' H74[l . In contrast to our FE-ED, both the hybridiza- 
tion parameters Vi and the site energies q are determined 
self-consistently in the Caffarel-Krauth scheme (CK-ED). 

In order to close the self-consistency cycle, a x 2 -fhting 
procedure |26I18| is used. Let T be a fictitious small tem- 
perature, typically T = 1/200, and ui n = (2n + l)nT are 
the corresponding Matsubara frequencies. For a given pa- 
rameter set (et, Vi) we solve the single-impurity model for 

the single-particle Green function Ga (w). Then, the self- 
consistency equation (|7§)l allows us to deduce the param- 
eter set for the next iteration from the minimization of 

2 / new T/now\ \ * Tl( n s) ft, , ■ ,new r/new\ 



n=0 



G^\iu n -e e ,V e ) . (123) 



Here um is a (large) upper cut-off. 




Fig. 16. Density of states for U = 6 as obtained from the 
l/U expansion (solid line), in comparison with two converged 
solutions of the Caffarel-Krauth Exact Diagonalization scheme 
(CK-ED) for n s = 11 sites (solid points). 
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In Fig. 1161 we show two typical converged results for 
the density of states for U = 6, starting from two different 
initial values for the set (Ve,ee). Note that the solutions 
are not unique, i.e., there are many self-consistent solu- 
tions for the insulator. At first glance, the two solutions 
appear to be very similar. The main peaks have almost 
the same weights and the same positions eg. This can 
be expected because the main peaks dominate x 2 m EHt • 
Moreover, the overall distribution of weights and energies 
by and large follows the true density of states from the 
analytical calculation. 

A closer look into the figures reveals a fundamental 
problem. In Fig. 116b , the first non- vanishing peak appears 
at Z\i n) (6)/2 = 0.966 whereas in Fig. \T§p the value for 

half the gap is estimated as Z\[ n) (6)/2 = 1.06. The same 
results follow if we use the difference of the ground-state 
energies for the definition of the gap. Obviously, we 
cannot give a unique value for /\ CK - ED ;(™s) 5 an d an ex- 
trapolation n s — > oo is impossible. The reason for this 
problem lies in the fact that the peak at the onset of the 
density of states is rather small so that its precise position 
and weight do not matter much in the x 2_n t (ESU - Thus, 
the onset energy may easily fluctuate by 10% or more, as 
seen in Fig. EH Therefore, the CK-ED approach cannot 
provide a reliable estimate for the gap. 

The failure of the self-consistency can be traced back 
to the fact that, for general filling, the Green function 
G( n ">(Lj) contains only n s poles with substantial weight 
whereas i?( ns )(u) is characterized by (2n s ) parameters 
(e e ,Vf) in the CK-ED. Too much flexibility in the hy- 
bridization function leads to non-unique solutions of the 
self-consistency equations. Our Fixed-Energy Exact Diag- 
onalization cures this problem because only n s parameters 
V? are to be determined. 



8 Conclusions 

In this work, we have studied the insulating phase of the 
half-filled Hubbard model on a Bethe lattice with infinite 
connectivity (band- width W = At). As long as the strong- 
coupling perturbation theory converges, the ground state 
of the Mott-Hubbard insulator has a finite entropy den- 
sity, s = ln(2), and the density of states of the lower and 
upper Hubbard bands increases as a function of frequency 
with the edge exponent a — 1/2, see (|57f) . Furthermore, 
we find that the high-order corrections to the gap as a 
function of 1/U are fairly small. Thus, the gap opens at 
the critical interaction strength £/j? c = 4.40 ±0.09. Our ex- 
plicit results to second order in 1/U i(5B|) provide a bench- 
mark test for analytical and numerical approaches to the 
Mott-Hubbard insulator in infinite dimensions. 

We have used our results from the 1/U expansion to 
set up a new numerical scheme for the effective single- 
impurity Anderson model in the Dynamical Mean-Field 
Theory. In our Fixed-Energy Exact Diagonalization (FE- 
ED), we start an outer self-consistency cycle with the po- 
sition and width of the Hubbard bands from perturba- 
tion theory which we resolve equidistantly to accuracy 



9t/(n s — 1) for n s < 15. The extrapolation n s — > oo 
thereby becomes systematic. In an inner self-consistency 
cycle we determine the hybridization strengths for given 
n s using a dynamical Lanczos procedure. From this we de- 
termine the gap after extrapolation to the thermodynamic 
limit. We merge the extrapolated gap from the 'energy cri- 
terion' and the 'weight criterion' in order to obtain a new 
estimate for the onset of the Hubbard bands. In this way, 
we iterate to convergence the outer self-consistency cycle. 

Our FE-ED very well reproduces the gap, the density 
of states, and the exponent a = 1/2 from perturbation 
theory. In this way we confirm our estimate for the transi- 
tion to J7 C FE - ED = 4.43±0.05. As a last check, we favorably 
compare our results with those from the Random Disper- 
sion Approximation which is an independent approach to 
the limit of infinite dimensions. Thus, we are confident 
that the strong-coupling perturbation theory, our Fixed- 
Energy Exact Diagonalization scheme for the Dynamical 
Mean-Field Theory, and the Random Dispersion Approx- 
imation provide a consistent and accurate description of 
the Mott-Hubbard insulator. 

Other analytical and numerical techniques for the solu- 
tion of the Dynamical Mean-Field Theory meet our bench- 
mark test with limited success. The best analytical ap- 
proximation is the Local Moment Approach which is quan- 
titatively reliable for the density of states down to U = 
5.5; it reproduces the correct exponent a = 1/2 but over- 
estimates the critical interaction strength, !7 C LMA = 4.8. It- 
erated Perturbation Theory reproduces the correct thresh- 
old exponent a = 1/2 but it becomes quantitatively and 
qualitatively unreliable below U = 6 and underestimates 
the stability of the Mott-Hubbard insulator, Ul PT = 5.2. 

Of the two exact-diagonalization (ED) schemes for the 
effective single-impurity problem existing in the literature, 
the two-chain approach fails because it leads to an ill- 
conditioned inverse problem, namely, the recovery of the 
density of states from its moments. The Caffarel-Krauth 
exact diagonalization approach (CK-ED) in star geometry 
gives a reasonable guess for the shape of the density of 
states, but it fails to provide a reliable estimate for the 
gap. For fixed system size, the fitting procedure results in 
different solutions with noticeable differences for the onset 
of the density of states. Therefore, it is not possible to set 
up a sensible extrapolation scheme for n s — > oo. 

In this work we have presented a qualitatively and 
quantitatively reliable analysis of the Mott-Hubbard in- 
sulator. It will provide a solid basis for further studies of 
the Mott-Hubbard metal-insulator transition. 
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A Iterated Perturbation Theory (IPT) in 
strong coupling 

We here write the IPT self-energy, cf. I|105[l - p08)l . in a 
form amenable to numerical calculations and addition- 
ally find its strong coupling form. We start by rearrang- 
ing 1)1121) as follows 



g a {u) = <? ba «V) + 



\m \ 



In order to calculate QE l7 (u)) for u> > 0, it is thus neces- 
sary, see 1)1290 , to perform only the two convolutions 1)132)1 
and 1)133)1 . ^sS a (u>) for lu < follows immediately by sym- 
metry, whence ^tS^^uu) follows by Hilbert transform. 

The above is what we implement numerically. How- 
ever, when U — > oo eq. (TL29I) simplifies. In this limit itlToll 
gives | mo | — > 1; note that |mrj| — > is not a self-consistent 
possibility, since \mo\ > 2/3, cf. I)117|) . Specifically, 



lu + irj sgn(w) 



r ba„d/ \ , r band, n , l m ol/ 2 , \ m o\/ 2 

+^| mo |(5(o;-0-)-<5(a;-0+)) , (124) 

where £/ ba _£ d (aj) are the retarded/advanced components 
of C/ band (o;). Inserting (1 12411 into (|108fl . and performing 
a simple contour integration, yields for II(uo) 



\mo\ = l -Jj2 



U -> oo 



(134) 



Since |mo| is the weight of the pole in ^(w), l|112|l . and 

J_ l^Grriu))] dcu = 7r, we have 



3£ ban V) - 0(£7- 2 ) , U -f oo 



(135) 



Now, eqs. (TT32T) and ifT^I result in /i(w) - C(J7" 4 ) and 
7 2 H ~ C(C/ -2 ), whereby (Xjgj) reduces to 



Z7( W ) 



(|m |/2) 2 (|m |/2) : 



w — i?7 cj + irj 

+ \mo\ {G h a ™ d (u>) - <£T»} + h(u) , (125) 



\m \ 



3£° and (w) + 



\m \ 



lu + ir] sgn(w) 



where 



(136) 

with corrections 0(U 2 ). Using 1)112)) and (|107|l this may 
be rewritten as 



W = - ^r^^i)^^! - w) . (126) ra(w) = 



-3|m | 2 l/ 2 /4 £/ 2 K| 3 /2 

EV2-w' + G ff (w) CT/2-w' 1 J ' 

(137) 

where lu' = lu + U/2. In the region of the lower Hubbard 
band only (lu' ~ 0(1)) eq. (|137|) may be expanded to give 



Note that S77(w) has a pole contribution at u; = 0, specif- 
ically 

I7(w) = y |m | 2 <5H + 7T band H , (127) 

with 

77 band H = |m | K^V) - #£ d (w)} + J iH • ( 128 ) 

Insertion of [fT2"4*|) and [fl"2"5|) into ifTUfj|) yields for the self- 
energy 



= -|-u/ + 3G CT (u;) 



(138) 



C7 



[w' 2 - 6u/G CT (w) + 3G ct (cj) 2 ] + 0(ET 



lu + irj sgn(w) 



+ I 3 (u;) + ^C/ 2 {/f( W ) - , (129) 



If self-energy terms of O (U 1 ) and above are ignored, then 
G&(lu) in the region of the lower Hubbard band reduces 
to 

Ga{uu) = - X 'l , (139) 
lu' — 2G cr (w) 

and the single-particle density of states becomes 



where 



/aH = U 2 r -prg^iu; - Wl )7T band ( Wl ) (130) 
J-oo 2vri 



^lmb.^M - ^-^JA-lu'Z , 

47T 



(140) 



and 



i r 00 dwiQfJi^i) 



Wl — lu =p 177 



(131) 



i.e., the correct strong coupling spectrum of a semi-ellipse 
with the full noninteracting band-width. Retaining the 
0(C/ _1 ) terms in eq. 1)138)1 gives the 0(U~ r ) corrections 
to the position and shape of the LHB. The resultant spec- 
trum is a distorted semi-ellipse of the form 



^sli (lu) may be found using a convolution involving only 
3<? band (o;), and likewise ^^(w) may be found from a con- 
volution of %g% and (Lu) and Q7I band (cj). For lu > 0, 



LHB 



fllPT V 



- / d Wl 3£ band ( Wl )3£ ban V- Wl ) ,(132) 



v J °^ 

l-(aiPT/UW + Q(U- 2 ) 

47T 



. (141) 



9/2(0,) = *L rd Wl Q£ band ( Wl )3i7 ban V- Wl ) . 

7f JO 



(133) 



where aipx = 5/4 and the LHB has band edges at lu' = 
(<1wt/U) ± 2. The IPT result may be compared with 
the exact result to this order which also has the form of 
eq. 1141)1 but with a different coefficient of a = 1/2. 
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B Local Moment Approach (LMA) in strong 
coupling 

The LMA of Logan et al. is described in detail in Ref. |25| . 
where it is proved that the LMA recovers the strong cou- 
pling single-particle spectrum up to and including terms 
of O(U ). Unlike IPT, the LMA also recovers the strong 
coupling spectrum in the antiferromagnetic phase |25I43| . 
The success of the LMA in strong coupling is unsurpris- 
ing since the LMA is explicitly motivated by ideas of hole 
motion in a spin background. In this appendix we extend 
the previous analysis to recover the 0(C/ _1 ) corrections 
to the LMA self-energy and single-particle spectrum in 
the strong-coupling limit. This is done in three stages. 
First the necessary equations of the LMA are given. Sec- 
ondly, it is shown that only the coupling of hole motion to 
zero-frequency spin-flips survive, with the coupling of hole 
motion to higher energy particle-hole excitations only en- 
tering at (D(U~ 2 ). Finally, the resultant self-energy is ex- 
panded to C(J7 _1 ) and the corresponding density of states 
is obtained. 

The bare propagators within the LMA are unrestricted 
Hartree Fock (UHF) propagators. Lattice sites are taken 
to have up-spin or down-spin permanent local moments 
at random, and labeled 'A-type' or 'B-type' respectively. 
We focus solely on an A-type site (Green functions for 
B-type sites follow by symmetry) and suppress the 'A' 
labels in what follows. For the Z = oo Bcthe lattice the 
UHF equations are simply 

G£M = -G°_ a (-uj) = [w + ^U | M | - GV)]- 1 , (142) 

where a = ±1 for ] / { respectively, and the symmetric 
single-particle Green function obeys G°(u>) — [G®(uj) + 
G|(w)]/2. The local moment is derived from 



duj[D^(uj) - Dl(w)] 

> 

D° a (uj) = -isgnH3G°H . 

The RPA on-site particle-hole propagator 

77+- H 



77+- M = 



i-un+-{u) 



with 



77+-H = - J —Gl(n)G°(n cj) 



(143) 
(144) 

(145) 

(146) 



has a zero-frequency pole |25| (as well as higher-energy 
Stoner bands) reflecting the presence of zero-frequency 
spin-flip excitations in the system. In the full LMA single- 
particle Green function, 



G a (u)= u + -U\fi\-G(uj)-E a (uj) 



(147) 



hole motion and accompanying spin-flips are dynamically 
coupled via the self-energy 



S ] {uj)=-S 1 {-u) = U 2 



d(2 
2~7ri 



zr+-(fl)gi(w + fl). 



(148) 

The LMA equations must be solved self-consistently, since 
the host Green function reads 

g <T (uj) = [uj + ~U\[i\-G(uj)]- 1 (149) 

with G{uS) = [G T (w) + G x H]/2. 

To calculate the LMA self-energy to ©(t/- 1 ), we pro- 
ceed in two stages. First we sketch the proof that only the 
zero-frequency spin-flip pole in 77 ^ (to) need be retained 
and the self-energy simplifies to 



u 



, i dwi ^x(^i) , ,o, r - i 



7T LO — LO\ — lTj 



(150) 

+ D(U-' 2 ) . (151) 



Subsequently we analyze (lo) to obtain the spectrum to 

oiu- 1 ). 

It has been shown [231 that the RPA particle-hole prop- 
agator may be written 



n+-(w) 



Q 



+ n+-(uj) + n+-(uj) , (152) 



where 77 + -(oj) and 77 + -(w) are the Stoner-like contribu- 
tions whose imaginary parts are bands centered at ±U 
respectively. The full self-energy is thus 



27 T (w) = QU 2 g^{u;) + S^{u;) 



(153) 



where 



SZf (w 



-U 2 dflV i (Q)^n+- {Q - uj) for uj > 



[U 2 J^dttV^tySII^-in- uj) forw<0 

(154) 

To show that (uj) ~ 0(U 2 ), it is necessary to examine 

77 H (uj) and hence 77q (uj) and G° a (uj). By expanding 
G\(lo) in the region of the LHB it quickly follows that 



\n\ = i-i/u 2 + o(u- 3 ) 



(155) 



Then 377 + (uj) may be shown from l(146[l to have a mi- 
nority band centered around uj ps — U with spectral weight 
[(1 - |/i|)/2] 2 = l/(4[/ 4 ) + 0(U- 5 ), and a majority band 
centered at uj « U with spectral weight [(1 + \fi\)/2} 2 = 
1 - 1/(U 2 ) + 0(U~ 3 ). The bands of the RPA propagator 
which follow from 1145(1 . 



Q77+-(cj) = 



577q + -H 

\i-un+-(u;)\ 



(156) 



are also centered around uj — ±U. Since 11^ (uj) ~ 0(1) 
in the region of the upper band it follows that the spec- 
tral weight in S77 + ~(w) is 0(U~ 2 ). In the region of the 
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lower Hubbard band (uj w —U/2) Kramers-Kronig rela- 
tions result in rni£~(u) = 1/(2U) + <D{IJ- 2 ), whence 
377+-(w) = 4%II+-(uj)+0(U- 5 ) and the spectral weight 
of (u>) is 1/U 4 . The remainder of the spectral weight 

of 377 H (uj) is contained in the pole weight Q. Eqs. i|145|) . 
(jl46|l may be used to show that 



n + -(cj -> oo) = n+-(u 



oo) = -— , (157) 



whence Kramers-Kronig relations give the missing weight 
as Q = 1 — 0(C/~ 2 ). Finally, since (in common with 

G®(uj)) has minority and majority spectral weights that 
approach 1/(2U 2 ) and 1 - 1/(2U 2 ) respectively, eq. (I154|l 
shows the contributions to QS^(u>) are 0(U~ 2 ). A more 
detailed analysis shows that in the region of the LHB 
S^(uj) is pure real and 0(U~ 4 ), and thus neglecting S^(u>) 
barely changes the LMA results throughout the insulating 
phase. Nevertheless the analysis outlined here is sufficient 
to arrive at (|150|) . 

Substituti ng iJTSOj l together with \fi\ = 1 - 1/U 2 + 
0(U~ 3 ) into fT4T| yields in the region of the LHB 



/ G T H-tf 2 sr>) + o(Er a ) 



where we have defined 



(158) 



uj'=uj + U/2, (159) 
and we have also used the LHB expansion 



which is correctly obtained within the LMA. The Hilbert 
transform l|151|) then yields 



1 



U%(uj) = -G T M + — i3G T H 
4t/ Jlhb w - wi 

= \ G]{uj) + w G ' (uj) ~ h + ° {u ~ 2) 



(164) 



where we have used 2ttD^(uj) = y/A — uj' 2 + 0(U~ 4 ) to 
simplify the C(J7 _1 ) correction. Eq. i|164l) is only valid 
in the region of LHB spectral density. Substituting 1164JI 
into (|158(l gives the single-particle Green function in this 
region, 



G T (w) = 

'(165) 

The corresponding single-particle spectrum in strong cou- 
pling including 1/U corrections is 



D T (w) = \U-[u)'- 



Qlma 
U 



o(u- 



y l-(a^MA/UW + Q(U' 2 ) (16g) 
2ir 

with the LMA coefficient olma = 3/4 instead of a = 1/2 
from the strong-coupling expansion. 



G{u)= l -G ] {u:)-± J + 0{U- 2 ) . (160) 

All that remains is to calculate Q7(uj) to 0(U~ 3 ). We 
focus exclusively on the LHB. First we expand Qi(uj) in 
powers of 1/U, 



i i 

u ~ u 2 
1 

u 3 



2U 



G(uj) 



2U 



G(uj) 



+ 0(U~ 4 ) .(161) 



Since in the region of the LHB (uj) = QQi(ui) by 
definition, we have 

U 2 %G7(uj) = ZG(uj) - j-%{[uj' - G(uj)] 2 } + 0(U~ 2 ) 

= ~3G T M - ±^{u,'G^w)} + 0(U- 2 ) , 

(162) 

where we have used <|16(Jll and the exact form of the single- 
particle Green function in strong coupling, 

G](uj) =uj'G^(uj)-1 + 0(U- 1 ) , (163) 



References 

1. N.F. Mott, Metal-Insulator Transitions, 2nd edition (Tay- 
lor and Francis, London, 1990). 

2. F. Gebhard, The Mott Metal-Insulator Transition (Sprin- 
ger, Berlin, 1997). 

3. A. Georges, G. Kotliar, W. Krauth, and M.J. Rozenberg, 
Rev. Mod. Phys. 68, (1996) 13. 

4. M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 
70, (1998) 1039. 

5. E.H. Lieb and F.Y. Wu, Phys. Rev. Lett. 20, (1968) 1445. 

6. F. Gebhard and A.E. Ruckenstein, Phys. Rev. Lett. 68, 
(1992) 244; F. Gebhard, A. Girndt, and A.E. Ruckenstein, 
Phys. Rev. B 49, (1994) 10926. 

7. W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, (1989) 
324. 

8. E. Miiller-Hartmann, Z. Phys. B 74, 507 (1989). 

9. U. Brandt and C. Mielsch, Z. Phys. B 75, (1989) 365, 
ibid. 79, (1990) 295. 

10. M. Jarrell, Phys. Rev. Lett. 69, (1992) 168. 

11. R.M. Noack and F. Gebhard, Phys. Rev. Lett. 82, (1999) 
1915. 

12. J. Schlipf, M. Jarrell, P.G.J, van Dongen, N. Bliimer, 
S. Kehrein, T. Pruschke, and D. Vollhardt, Phys. Rev. 
Lett. 82, (1999) 4890. 

13. M.J. Rozenberg, R. Chitra, and G. Kotliar, Phys. Rev. 
Lett. 83, (1999) 3498. 

14. W. Krauth, Phys. Rev. B 62, (2000) 6860. 



22 



M.P. Eastwood et al.: Mott-Hubbard Insulator in Infinite Dimensions 



15. R. Bulla, Phys. Rev. Lett. 83, (1999) 136. 

16. R. Bulla, T.A. Costi, and D. Vollhardt, Phys. Rev. B 64, 
045103 (2001). 

17. R. Bulla and M. Potthoff, Eur. Phys. J. B 13, (2000) 257. 

18. Y. Ono, R. Bulla, A.C. Hewson, and M. Potthoff, Eur. 
Phys. J. B 22, (2001) 283. 

19. G. Moeller, Q. Si, G. Kotliar, M.J. Rozenberg, and D.S. 
Fisher, Phys. Rev. Lett. 74, (1995) 2082. 

20. D.E. Logan and P. Nozieres, Phil. Trans. R. Soc. London 
A 356, (1998) 249; P. Nozieres, Eur. Phys. J. B 6, (1998) 
447. 

21. S. Kehrein, Phys. Rev. Lett. 81, (1998) 3912. 

22. A. Georges and G. Kotliar, Phys. Rev. Lett. 84, (2000) 
3500; S. Kehrein, Phys. Rev. Lett. 84, (2000) 3501. 

23. M.J. Rozenberg, G. Kotliar, and X.Y. Zhang, Phys. Rev. 
B 49, (1994) 10181. 

24. T. Pruschke, D.L. Cox, and M. Jarrell, Europhys. Lett. 21, 
(1993) 593; Phys. Rev. B 47, (1993) 3553; M. Jarrell, 
J.K. Freericks, and T. Pruschke, Phys. Rev. B 51, (1995) 
11704. 

25. D.E. Logan, M.P. Eastwood, and M.A. Tusch, J. Phys. 
Cond. Matt. 9 (1997) 4211. 

26. M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, (1994) 
1545. 

27. M.J. Rozenberg, G. Moller, and G. Kotliar, Mod. Phys. 
Lett. B 8 (1994), 535. 

28. See also Q. Si, M.J. Rozenberg, G. Kotliar and A.E. Ruck- 
enstein, Phys. Rev. Lett. 72, (1994) 2761; M.J. Rozen- 
berg, G. Kotliar, H. Kajiiter, G.A. Thomas, D.H. Rapkine, 
J.M. Honig, and P. Metcalf, Phys. Rev. Lett. 75, (1995) 
105. 

29. E. Kalinowski and F. Gebhard, J. Low Temp. Phys. 126, 
(2002) 979. 

30. P.G.J, van Dongen, F. Gebhard, and D. Vollhardt, Z. Phys. 
B 76, (1989) 199. 

31. E. Economou, Green's Functions in Quantum Physics, 2nd 
edition (Springer, Berlin, 1983). 

32. J. Hubbard, Proc. Roy. Soc. London Ser. A 276, (1963) 
238; ibid. 277, (1963) 237. 

33. A.L. Fetter and J.D. Walecka, Quantum Theory of Many- 
Particle Systems (McGraw-Hill, New York, 1971). 

34. T. Kato, Prog. Theor. Phys. 4, (1949) 154; A. Messiah, 
Quantum Mechanics, Vol. 2 (North-Holland, Amsterdam, 
1962) § 16. 

35. M. Takahashi, J. Phys. C 10, (1977) 1289. 

36. P.W. Anderson, Phys. Rev. 115, (1959) 2. 

37. W. Metzner, P. Schmit, and D. Vollhardt, Phys. Rev. B 45, 
(1992) 2237; W.F. Brinkman and T.M. Rice, Phys. Rev. 
B 2, (1970) 1324. 

38. P.G.J, van Dongen and D. Vollhardt, Phys. Rev. Lett. 65, 
(1990) 1663; P.G.J, van Dongen, Phys. Rev. B 45, (1992) 
2267. 

39. M. Potthoff, preprint cond-mat/0301137 (2003), unpub- 
lished. 

40. A. Georges and G. Kotliar, Phys. Rev. B 45, (1992) 6479. 

41. M.P. Eastwood, Ph.D. thesis, Oxford University (1998), 
unpublished. 

42. J. Hubbard, Proc. Roy. Soc. London 281, (1964) 401. 

43. D.E. Logan, M.P. Eastwood, and M.A. Tusch, Phys. Rev. 
Lett. 76, (1997) 4785. 



